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Abstract. Smale's a-theory uses estimates related to the convergence of Newton's 
method to certify that Newton iterations will converge quadratically to solutions to a 
square polynomial system. The program alphaCertified implements algorithms based 
on a-theory to certify solutions of polynomial systems using both exact rational arith- 
metic and arbitrary precision floating point arithmetic. It also implements algorithms 
that certify whether a given point corresponds to a real solution, and algorithms to 
heuristically validate solutions to overdetermined systems. Examples are presented to 
demonstrate the algorithms. 



Introduction 

Current implementations of numerical homotopy algorithms [TJ I32J [38] such as PHC- 
pack jH], HOM4PS [27], Bertini [2], and NAG4M2 [28] routinely and reliably solve sys- 
tems of polynomial equations with dozens of variables having thousands of solutions. Here, 
'solve' means 'compute numerical approximations to solutions.' In each of these software 
packages, the solutions are validated heuristically — often by monitoring iterations of New- 
ton's method. This works well in practice, giving solutions that are acceptable in most 
applications. However, a well-known shortcoming of numerical methods for computing 
approximate solutions to systems of polynomials is that the output is not certified. This 
restricts their use in some applications, including those in pure mathematics. The program 
alphaCertified is intended to remedy this shortcoming. 

In the 1980 's, Smale [5B] and others investigated the convergence of Newton's method, 
developing a-theory Ch. 8]. This refers to a computable positive constant a(f,x) 
depending upon a system / : C n — > C n of polynomials and a point x G C n such that, if 

a(f,x) < 13 ~ 4 3v/I ^ « 0.157671, 

then iterations of Newton's method starting at x will converge quadratically to a solution 
to /, which is a point £ G C n with /(£) = 0. In principle, Smale's a-theory provides 
certificates for validating numerical computations with polynomials. 
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Current implementations of numerical homotopy algorithms do not incorporate a- 
theory to certify their output or their path-tracking. There have been two projects which 
use fixed double precision and focus on certified path-tracking. Malajovich [30] released 
the most recent version of his Polynomial System Solver in 2003, which uses a-theory 
to certify toric path-tracking algorithms, but he states that "[it] is actually not intended 
for an end user." Beltran and Leykin [5] have recently shown how to use a-theory to 
certify path-tracking, and hence the output of numerical homotopy algorithms. While 
they demonstrate that certification can dramatically affect the speed of computation, this 
is an important development, as certified path-tracking is necessary for applications such 
as numerical irreducible decomposition [37] or computing Galois groups [29]. They are 
continuing this line of research. 

We describe a program, alphaCertified, that implements elements of a-theory to 
certify numerical solutions to systems of polynomial equations using both exact rational 
and arbitrary precision floating point arithmetic. As it only certifies the output of a 
numerical computation, it avoids the bottlenecks of certified tracking, while delivering 
some of its benefits. Given a square polynomial system / : C n — > C", alphaCertified uses 
Smale's a-theory to answer the following three questions for a finite set of points X C C n : 

(1) From which points of X will Newton's method converge quadratically to some 
solution to /? 

(2) From which points of X will Newton's method converge quadratically to distinct 
solutions to /? 

(3) If / is real ({fi, . . . , f n } = {/i, . . . , f n }), from which points of X will Newton's 
method converge quadratically to real solutions to /? 

Often, a sharp upper bound B on the number of roots to a square polynomial system / 
is known. Given a set of B points, alphaCertified can be used to certify that iterations 
of Newton's method starting from each point in the set converge quadratically to some 
solution to / and that these solutions are distinct. This guarantees that each of the B 
roots of / can be approximated to arbitrary accuracy using Newton's method. Moreover, 
alphaCertified can certify how many of the B solutions to / are real when / is real. 

A polynomial system / : C n — > C N is overdetermined ii N > n, that is, if the number 
of polynomials exceeds the number of variables. Dedieu and Shub [12] studied Newton's 
method for overdetermined polynomial systems and gave conditions which guarantee qua- 
dratic convergence of its iterations. Unlike square systems, the fixed points of this overde- 
termined Newton's method need not be solutions. For example, x = 1 is a fixed point of 

x 

Newton's method applied to f(x) = ^ . 

The program alphaCertified validates solutions to overdetermined systems. Given a 
finite set IcC" and an overdetermined system, it generates two or more random square 
subsystems, answers the three questions above for each, and compares the results. In 
particular, given 5 > 0, it can certify that, for a given approximate solution to two or 
more random subsystems, the associated solutions all lie within a distance 5 of each other. 
For a given S, this heuristically validates solutions to overdetermined systems. 
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In summary, alphaCertified is novel in each of the following ways. It implements 
algorithms from a-theory using either exact rational or arbitrary precision floating point 
arithmetic. When using exact rational arithmetic with a square polynomial system, its 
implementation of a-theory is completely rigorous. It certifiably determines if an ap- 
proximate solution corresponds to a real solution, which may be used to count the real 
solutions to a polynomial system, and it uses a-theory to obtain information on the roots 
of overdetermined systems. The examples we give demonstrate the practicality of certifi- 
cation based on a-theory, and its viability as an alternative to exact symbolic methods, as 
the certificates for square systems when using exact rational arithmetic are mathematical 
proofs of computed results. 

In Section [TJ we review the concepts of a-theory utilized by alphaCertified. Section [2] 
presents the algorithms for square polynomial systems while Section [3] describes our ap- 
proach to overdetermined polynomial systems. Implementation details are presented in 
Section H] with examples presented in Section verifying some computational results in 
kinematics and generating evidence for conjectures in enumerative real algebraic geometry. 



1. SMALE'S a-THEORY 

We summarize key points of Smale's a-theory for square polynomial systems that are 
utilized by alphaCertified. More details may be found in 0, Ch. 8]. 

Let / : C n — > C n be a system of n polynomials in n variables with common zeroes 
V(/) := {£ E C n | /(£) = 0}, and let Df(x) be the Jacobian matrix of the system / at x. 
Consider the map Nf : C n — > C n defined by 



N f {x) 



x — Df(x) 1 f(x) if Df(x) is invertible, 
x otherwise. 



The point Nf(x) is called the Newton iteration of f starting at x. For k EN, let 

N'f(x) := N f o ■ • • o N f (x) 
k times 

be the k th Newton iteration of / starting at x. 

Definition 1. Let / : C" — > C™ be a polynomial system. A point x E C n is an approximate 
solution to / with associated solution £ E V(f) if, for every k E N, 

(i) \\N h f {x)-a < Q) 2 \\x-z\\. 

That is, the sequence {Nj(x) \ k E N} converges quadratically to £. Here, || ■ || is the usual 
hermitian norm on C n , namely . . . ,x n )\\ = (\xi\ 2 + ■ ■ ■ + la^l 2 ) 1 / 2 . 

Smale's a-theory describes conditions that imply a given point x is an approximate 
solution to /. It is based on constants a(f, x), (3(f, x), and j(f, x). If Df(x) is invertible, 
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these are 



(2) 7(/» : = sup 

k>2 



<x{f,x) := /3{f,x)^(f,x) , 

/3(f,x) := \\x-N f (x)\\ = \\Df(x)- l f(x)\\, and 
Df{xY l D k }\x) 



k\ 



If x G V(/) is such that Df(x) is not invertible, then we define a(f,x) := /3(f,x) := 
and j(f,x) := oo. Otherwise, if x ^ V(/) and Df(x) is not invertible, then we define 
«(/, a:) := /?(/, := 7(/» x ) := °°- 

In the formula ([2]) for j(f,x), the derivative D k f(x) [251 Chap. 5] to / is the 
symmetric tensor whose components are the partial derivatives of / of order k. It is a 
linear map from the fc-fold symmetric power S k C n of C n to C n . The norm in ([2]) is the 
operator norm of D f (x) _1 D k f '(%) : S k C n —> C n , defined with respect to the norm on S k C n 
that is dual to the standard unitarily invariant norm on homogeneous polynomials [25] , 



£kl7C), 



where v — (yi, . . . , v n ) is an exponent vector of non-negative integers with X — X ^ * X yP^ - 
= i^i + • • • + v n , and (J = i is the multinomial coefficient. 

The following version of Theorem 2 from page 160 of [9] provides a certificate that a 
point x is an approximate solution to /. 

Theorem 2. If f : C n — >• C n is a polynomial system and x G C n u>z£/j 



1 Q _ Q,/1 7 

(3) «(/,*) < ^— « 0.157671, 

then x is an approximate solution to f. Additionally, \\x — £|| < 2(3(f,x) where £ G V(/) 
zs t/ie associated solution to x. 

Remark 3. If a(f, x) > |, then x may not be an approximate solution to /. For example, 
for f(x) = x 2 , if x 7^ 0, then x is not an approximate solution to / yet a(f, x) = \. 

For a polynomial system / : C n — > C™ and a point x G C n , we say that x is a certified 
approximate solution to / if (j3J) holds. 

Theorem 4 and Remark 6 of [HI Ch. 8] give a version of Theorem [2] that alphaCerti- 
fied uses to certify that two approximate solutions have the same associated solution. 

Theorem 4. Let f : C n — > C n be a polynomial system, x G C n with a(f,x) < 0.03 and 
£ G V(/) f/ie associated solution to x. If y G C n wift 

1 

||x — j/ < 



20 7 (/,x)' 

t/ien j/ 2s an approximate solution to f with associated solution £. 
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1.1. Bounding higher order derivatives. The constant j(f,x) encoding the behavior 
of the higher order derivatives of / at x is difficult to compute, but it can be bounded 
above. For a polynomial g : C n — > C of degree d, say g = 'Yl\ u \ <d a u x u , define 

2 x!(<f-|"l)! 



d\ 



\v\<d 

Then || ■ || is the standard unitarily invariant norm on the homogenization of g. For a 
polynomial system / : C™ — > C n , define 



^H/ill 2 where f(x) 



i=l 



fn(x) 



and for a point ieC", define 



n 



\x\\ 2 = l + 



i=l 



Let A^)(x) be the n x n diagonal matrix with 

A (d) (x) M := d^W?" 1 , 
where is the degree of /j. If Df{x) is invertible, define 

M/,x) := max{l, ll/H ■|| J D/(x)- 1 A (d) (x)||}. 
The following version of Proposition 3 from §1-3 of [35] gives an upper bound for j(f, x). 

Proposition 5. Let f : C n — > C n be a polynomial system with d{ = deg fi and D = max d{. 
If x G C n such that Df(x) is invertible, then 

(4) 7 (/.«) < 



2||x||i 

2. Algorithms for square polynomial systems 

Let / : C™ — » C™ be a square polynomial system and X = {x\, . . . , Xk} C C™ be a 
set of points. We describe the algorithms implemented in alphaCertified which answer 
the three questions posed in the Introduction. These algorithms are stated for a polyno- 
mial system with complex coefficients, but are implemented for polynomial systems with 
coefficients in Q[v~T] using both exact and arbitrary precision arithmetic. 

For each i = l,...,k, alphaCertified first checks if /(xj) = 0. If f(xi) ^ 0, then 
alphaCertified determines if Df{xj) is invertible. If it is, alphaCertified computes 
f3(f,Xi) and upper bounds for a(f,Xi) and y(f,Xi) using the following algorithm. 

Procedure (a,j3, 7) = ComputeConstants(/, x): 

Input: A square polynomial system / : C n — > C n and a point x € C n such that 
Df(x) is invertible. 

Output: a := /3 ■ 7, j3 := \\Df(x)~ 1 f(x) ||, and 7, where 7 is the upper bound for 
y(f,x) given in Proposition [5] 
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The next algorithm uses Theorem [2] to compute a subset Y of X containing points that 
are certified approximate solutions to /. 
Procedure Y = CertifySolns(/, X): 

Input: A square polynomial system / : C n — > C n and a set X = {xi, . . . , Xk} C C n . 

Output: A set Y C X of approximate solutions to /. 

Begin: 

(1) Initialize Y := {}. 

(2) For j — 1, 2, . . . , k, if f(xj) = 0, set Y := YU{xj}, otherwise, do the following 
if Df(xj) is invertible: 

(a) Set (a,/3, 7) := ComputeConstants(/, x^-). 

(b) If a < 13 ~ 4 3 ^ , set y := F U {x,}. 
Return: y 

As alphaCertified uses the upper bound for j(f,x) of Proposition [51 it may fail to 
certify a legitimate approximate solution x to /. In that case, a user may consider retrying 
after applying a few Newton iterations to x. The software alphaCertified does not 
invoke an automatic refinement to inputs that it does not certify. This is because Newton 
iterations may have unpredictable behavior (attracting cycles, chaos) when applied to 
points that are not in a basin of attraction. However, alphaCertified does provide the 
functionality for the user to do this refinement. 

Suppose that x is an approximate solution to / with associated solution £ such that 
Df(£) is invertible. Since x is an approximate solution, /?(/, iVj?(x)) converges to zero. 
Since j(f, x) is the supremum of a finite number of continuous functions of x, j(f, Nj(x)) 
is bounded. In particular, a(f,Nf(x)) converges to zero. 

Given approximate solutions X\ and x<i to / with associated solutions £1 and £2, respec- 
tively, Theorems [2] and H] can be used to determine if £1 and £2 are equal. In particular, if 

||xi-x 2 || > 2(/3(/,x 1 ) + /3(/,x 2 )), 
then £1 7^ £ 2 by Theorem [21 If on the other hand we have 

a(f,Xi) < 0.03 and ||xi-x 2 || < on /, r 

2(W, Xi) 

for either i = 1 or i = 2, then £1 = £2 by Theorem HI This justifies the following algorithm 
which determines if two approximate solutions correspond to distinct associated solutions. 
Procedure isDistinct = CertifyDistinctSoln(/, x\, x 2 ): 

Input: A square polynomial system / : C n — > C n and approximate solutions X\ and 
x 2 to / with associated solutions £1 and £2, respectively, such that Df(£i) and 
Df(^2) are invertible. 

Output: A boolean isDistinct that describes if £1 7^ £2- 

Begin: Do the following: 

(a) For % = 1,2, set (ci:i,/3i, 74) := ComputeConstants(/, xi). 

(b) If ||xi - 3 2 1| > 2(/?i + /3 2 ), Return True. 

(c) If a, < 0.03 and ||xi — x 2 || < , for either i = 1 or i = 2, Return False. 

w 11 1 iii 207 . ; 
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(d) For i = 1,2, update Xi := Nf(xi) and return to (a). 
This will halt, determining whether or not £1 = £2 as /?(/, Nf(xi)) decreases quadratically 
with k, while ^(f, Nf(xi)) is bounded. 

2.1. Certifying real solutions. A polynomial system / : C n — y C n is realii . . . , / n } = 
{A) • • • 5 /n}- in that case, solutions to f(x) = are either real or occur in conjugate pairs. 
Also, Nf(x) = Nf(x) for x G C n so that Ay: M n — > W 1 is a real map. Theorems [2] and H] 
can be used to determine if an approximate solution of / is associated to a real solution. 
Let x be an approximate solution to / with associated solution £. We do not assume 
that x is real, for numerical continuation solvers yield complex approximate solutions. By 
assumption, x is also an approximate solution to / with associated solution £. If 

> 2 (/?(/, x) + 0(f,x)) = W,x), 

then £ 7^ £ by Theorem [2] since 

||£-£|| > ||z-z||- > 0. 
Consider the natural projection map 7r R : C n — )• M n defined by 

/ \ x + x 

M x ) = 

Since ||x — x\\ = 2\\x — 7T]r(x)||, £ is not real if 

(5) ||x-7r R (x)|| > 2/3(f,x). 

We have both a local and a global approach to show that £ is real. For the local ap- 
proach, Theorem H] implies that 7Tm(x) is also an approximate solution to / with associated 
solution £ if 

(6) a(f,x) < 0.03 and \\x - n R (x)\\ < ) . 

207(/, x) 

Since Nj is a real map and ir^(x) G M n , this implies that £ G M. n . 

We could also have showed that both x and x correspond to the same solution to deduce 
that £ = £. If 

a(f,x) < 0.03 and \\x — x\\ < 77 — r, 

207(J, x) 

then Theorem H] implies that £ = £. This is more restrictive then ([6]) since ||x — x\\ = 

2\\X - 7T R (x)||. 

When a(f,x) < 0.03, (JSj) and (J6]) yield closely related statements. Since 
h Mf ^ 5a(/,x) 5-0.03 1 



3" w ' ' 3 7 (/,x) 3 7 (/,x) 20 7 (/,x)' 

we know that £ is real if \\x — 7r R (x)|| < |/3(/, x) and not real if ||x — vr K (x)|| > 2f3(f, x). 

The following algorithm uses the local approach of (jSJ) and (EJ) to determine if an 
approximate solution corresponds to a real associated solution. 

Procedure isReal = Certify RealSoln(/, x): 
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Input: A real square polynomial system / : C n — > C n and an approximate solution 

x G C n with associated solution £ such that -D/(£) is invertible. 
Output: A boolean isReal that describes if £ G W 1 . 
Begin: Do the following: 

(a) Set (a,/3, 7) := ComputeConstants(/, x). 

(b) If \\x — 7r R (x)|| > 2/3, Return False. 

(c) If a < 0.03 and \\x — 7Tr(2;)|| < , Return True. 

(d) Update x := Nf(x), and return to (a). 

For the global approach to certifying real solutions, suppose that we know a priori 
that / has exactly k solutions. Suppose that X\ j . . . , Xfe £1X6 approximate solutions of / 
with distinct associated solutions. If, for all j ^ i, xi and Xj also correspond to distinct 
solutions, then xi and xl must correspond to the same solution, which is therefore real. 
This global approach requires a priori knowledge about V(f) as well as approximate 
solutions corresponding to each solution to /. While it cannot be applied to all systems, 
it is an alternative to the test based on j(f, x). 



2.2. Certification algorithm. For a given set of points X and a polynomial system /, 
CertifySolns, CertifyDistinctSoln, and CertifyRealSoln answer the three questions 
posed in the Introduction. We provide a sketch of the algorithm. 

Procedure (A,D,R) = CertifyCount(/, X): 

Input: A square polynomial system / : C n — > C" and a finite set of points X = 
{xi, . . . , xe} C C n such that if Xj is an approximate solution with associated solu- 
tion then Df(£j) is invertible. 

Output: A set A C X consisting of certified approximate solutions to /, a set 
D C A consisting of points which have distinct associated solutions, and, if / is a 
real map, a subset R C D consisting of points which have real associated solutions. 

Begin: 

(1) Set A := CertifySolns (f,X). 

(2) Set riA '■= \A\ and enumerate the points in A as ai, . . . , a nA . 

(3) For j = 1, . . . , tia, set Sj := True. 

(4) For j = 1,...,tia and for k = j + 1,...,ua, if Sj and Sk are True, set 
Sk '■= CertifyDistinctSoln(/, cij, a&). 

(5) Set D := {aj \ Sj = True}. 

(6) Initialize R := {}. 

(7) If / is a real polynomial system, do the following: 

(a) Set nr> '■= \D\ and enumerate the points in D as d\, . . . , d nD . 

(b) For j = l,...,n D , if CertifyRealSoln(/, dj) is True, update R : = 
RU{dj}. 
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3. OVERDETERMINED POLYNOMIAL SYSTEMS 



When N > n, the polynomial system / : C n — > C is overdetermined. Dedieu and 
Shub [12j studied the overdetermined Newton's method whose iterates are defined by 



where Df'(x) is the Moore-Penrose pseudoinverse of Df(x) [T71 § 5.5.4] to determine 
conditions that guarantee quadratic convergence. Since the fixed points of Nf may not 
be solutions to the overdetermined polynomial system /, this approach cannot certify 
solutions to overdetermined polynomial systems. 

We instead certify that points are associated solutions to two or more random square 
subsystems using the algorithms of Section [2j An additional level of security may be 
added by certifying that, for a given point which is an approximate solution to two or 
more random square subsystems, the associated solutions lie within a given distance of 
each other. As with the overdetermined Newton's method (JTj), this also cannot certify 
solutions to overdetermined polynomial systems, which is still an open problem. 

Let R: C N — > C™ be a linear map, considered as a matrix in C nxN . Then 1Z(f)(x) = 
Ro f(x) gives a square polynomial system 71(f) : C n — > C n . Since V(f) C V (71(f)) for 
any R, we call 71(f) a square subsystem of /. There is a nonempty Zariski open subset 
A C C nxN such that for every R G A and every x G V(f), null Df(x) = {0} if and only if 
D7Z(f)(x) is invertible. Moreover, for every x G V (71(f)) \ V(f), D7Z(f)(x) is invertible. 
See [38] for more on square subsystems 71(f). 

Define L = {f(x) \ x G C™} C which has dimension at most n possibly passing 
through the origin. A dimension-counting argument yields that there is a nonemtpy 
Zariski open set B C A x A C C nxN x C nxN such that, for every (R U R 2 ) G B, K = 
null R\ n null R 2 C is a linear space of dimension max{iV — 2n, 0} passing through 
the origin and K H L C {0}. In particular, if IZi(f) — Ri o /, then 



In addition, suppose that x is an approximate solution to both 7t\(f) and 7Z 2 (f) with 
associated solutions £i and respectively. For k G N, define = N^^(x) for z = 1,2. 
If £i 7^ there exists G N such that 



certifying that ||^i — < 5- In particular, this certifies that the solutions £i and ^2 to 
TZi and 7Z 2 associated to the common approximate solution x lie within a distance 5 of 
each other. For (5< 1, this heuristically shows that £1 = £2- 

In summary, if x is a certified approximate solution to two different square subsystems 
with distinct associated solutions, a certificate can be produced demonstrating this fact. 
Also (but not conversely), for any given tolerance 5 > 0, a certificate can be produced 



(7) 



N f (x) := x-Df(xYf(x), 



V(7Z 1 (f))nV(7Z 2 (f)) 



V(f). 
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that the distance between the associated solutions to the two square subsystems is smaller 
than 5. 

An additional test using the function residual could be added to this process. The 
following lemma describes such a test. 

Lemma 6. Let f : C n — > C N be an overdetermined polynomial system, R G C nxN , 
and x be an approximate solution to IZ(f) :— Ro f with associated solution £ such that 
a(TZ(f),x) < 0.0125. Then there exists e > such that if there exists y G C n satisfying 

\\x — v\\ < ; „. — r and \\f(v)\\ < e, 

thence V(f). 

Proof. Define v = — 7^7j\ — \ anc ^ v ) = {v ^ ^ N I ll x ~~ v\\ — u }- ^ e n °te that if 

j(lZ(f),x) = oo, since B(x, v) = {x}, it is easy to verify that we can take 

e = f 1 if f(x) = 0, 
\ ii^^H otherwise. 

Hence, we can assume that 7(7?.(/),x) < oo. Since 

lb- £11 < 2B(K(f) x) = 2 ^m^l < °- 025 = v 

£ G B(x, v). Moreover, Theorem H yields that B(x, v) n V{K(f)) = {£}• 

Assume £ ^ V(/). Since V(f) C V(TZ(f)), B(x,v) n V(/) = 0. In particular, = 
||/(^)|| is positive on the compact set B(x, u). Thus, there exists e > such that ||/(|/)|| > 
e for all y G B(x, v). □ 

Remark 7. For Lemma M to give an algorithm, we would need a general bound for the 
minimum of a positive polynomial on a disk. In cases when such a bound is known, e.g., 
j, it is too small to be practical. 

4. Implementation details for alphaCertified 

The program alphaCertified is written in C and depends upon GMP [19] and MPFR 
libraries to perform exact rational and arbitrary precision floating point arithmetic. 
When using rational arithmetic, all internal computations are certifiable. Because of the 
bit length growth of rational numbers under algebraic computations, alphaCertified al- 
lows the user to select a precision and use floating point arithmetic in that precision to 
facilitate computations. Since floating point errors from internal computations are not 
fully controlled, alphaCertified only yields a soft certificate when using the floating point 
arithmetic option. When the polynomial system is overdetermined, alphaCertified dis- 
plays a message informing the user about what it has actually computed. 

Three input files are needed to run alphaCertified. These files contain the polynomial 
system, the list of points to test, and the user-defined settings. See [22J for more details 
regarding exact syntax of these files. The polynomial system is assumed to have rational 
complex coefficients and described in the input file with respect to the basis of monomials. 
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That is, the user inputs the coefficient and the exponent of each variable for each monomial 
term in each polynomial of the polynomial system. 

The set of points to test are assumed to have either rational coordinates if using rational 
arithmetic or floating point coordinates if using floating point arithmetic. When using 
floating point arithmetic, the points are inputted in the precision selected by the user. 

The list of user-defined settings includes the choice between rational and floating point 
arithmetic, the floating point precision to use for the basic computations if using floating 
point arithmetic, and which certification algorithm to run. The user can also define 
a value, say r > 0, such that, for each certified approximate solution, the associated 
solution will be approximated to within 10~ r and printed to a file. 

The specific output of alphaCertified depends upon the user-defined settings. In each 
case, an on-screen table summarizes the output as well as a file that contains a human- 
readable summary for each point. The other files created are machine- readable files that 
can be used in additional computations. 

Linear solving operations are performed using an LU decomposition and the spectral 
matrix norm is bounded above using the Frobenius norm. This choice further worsens the 
approximation of 7 described in Proposition El which has two direct consequences on the 
performance of the algorithms. First, this requires that the value of (3 must be smaller 
in order to certify approximate solutions in CertifySolns. Second, algorithms Certify- 
DistinctSoln and Certify RealSoln may need to utilize extra Newton iterations. Apart 
from the added computational cost, the use of GMP and MPFR allows alphaCertified to 
still perform these computations even when using such an approximation of 7. 

When using rational arithmetic, alphaCertified avoids taking square roots when test- 
ing the required inequalities. When using floating point arithmetic, as an effort to con- 
trol the floating point errors, the internal working precision is increased when updating 
the point via a Newton iteration, for instance in Step (d) of CertifyDistinctSoln and 
Step (d) of CertifyRealSoln. 

The software alphaCertified determines if a square polynomial system / in n variables 
is real using two tests. The first test determines if the coefficients of / are real. The 
second selects a pseudo-random point y G Q n and determines if {fi(y), ■ ■ ■ , fn(y)} = 

{h(y),---,fn(y)}- 

The user either instructs alphaCertified to bypass all tests and declare that / is real, 
or which tests to use. If all tests fail, then alphaCertified bypasses the real certification. 
Otherwise, for each approximate solution x with associated solution £, alphaCertified de- 
termines if there exists a real approximate solution that also corresponds to £. If the user 
incorrectly identified / as real, then £ may not be real. Therefore, alphaCertified dis- 
plays a message informing the user about what it actually has certified. 

For an overdetermined polynomial system /, alphaCertified only checks to see if all 
of the coefficients of / are real. In this case, alphaCertified randomizes / using real 
matrices to obtain real square subsystems. 
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5. Computational examples 

We used alphaCertified to study four polynomial systems whose number of real solu- 
tions is relevant. Two are from kinematics and two are from enumerative geometry. All 
involve polynomial systems that are not easily solved using certified methods from sym- 
bolic computation. The files used in the computations, as well as instructions for their 
use, are found on our website [22] • Computations of Sections 15 .31 and I5T41 used nodes of the 
Brazos cluster [TU] that consist of two 2.5 GHz Intel Xeon E5420 quad-core processors. 

5.1. Stewart-Gough platform. The Stewart-Gough platform is a parallel manipulator 
in which six variable-length actuators are attached between a fixed frame (the ground) 
and a moving frame (the platform) [HI 00]. Each position of the platform uniquely 
determines the lengths of the six actuators. However, the lengths of the actuators do not 
uniquely determine the position and orientation of the platform, as there are typically 
several assembly modes, called positions. 

A generic platform with generic actuator lengths has 40 complex assembly modes. 
Dietmaier [13] used a continuation method to find a platform and leg lengths for which 
all 40 positions are real. While his formulation as a system of polynomial equations 
and conclusions about their solutions being real have been reproduced numerically (this 
is a problem in Verschelde's test suite [12]), these computations only give a heuristic 
verification of Dietmaier's result. 

We modified the polynomial system from Verschelde's test suite, which uses the parame- 
ters obtained by Dietmaier, by converting the floating point numbers to rational numbers. 
We then ran PHCpack [41] on the resulting polynomial system to obtain 40 numerical 
solutions to the system, each of which it identified as real. After converting the floating 
point coordinates of the solutions to rational numbers, we ran alphaCertified using these 
rational polynomials and rational points. It verified that these 40 points correspond to 
distinct real solutions. This gives a rigorous mathematical proof of Dietmaier's result. 

5.2. Four-bar linkages. A four-bar linkage is a planar linkage consisting of a triangle 
with two of its vertices connected to two bars, whose other endpoints are fixed in the 
plane. The base of the triangle, the two attached bars, and the implied bar between the 
two fixed points are the four bars. 



A general linkage has a one-dimensional constrained motion during which the joints may 
rotate, and the curve traced by the apex of the triangle is its workspace curve. 

The nine-point path synthesis problem asks for the four-bar linkages whose workspace 
curve contains nine given points. Morgan, Sommese, and Wampler [33J used homotopy 
continuation to solve a polynomial system which describes the four-bar linkages whose 
workspace curves pass through nine given points. They found that for nine points V = 
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{P , . . . , P 8 } C C 2 in general position, there are 8652 isolated solutions. Due to a two- 
fold symmetry, there are 4326 distinct four-bar linkages which appear in 1442 triplets, 
called Roberts cognates. We used alphaCertified to produce a soft certificate that the 
polynomial system has at least 8652 isolated solutions and, for a specific set of nine real 
points, certified the number of real solutions among these 8652 solutions. 

If V C M 2 , the formulation of [33] is not a real polynomial system. The usual approach 
of writing the variables using real and imaginary parts gives a real polynomial system 
f-p consisting of four quadratic and eight quartic polynomials. For nine points V = 
{Po, . . . , Ps} C C 2 , the polynomial system f-p : C 12 —> C 12 depends upon the variables 

{ai, a 2 , m, n 2 , xi, x 2 , h, b 2 , mi, m 2 , yi, y 2 } ■ 

Define the complex numbers 

a = ax + v^— T • a 2 , n = ri\ + v^— T • n 2 , x = X\ + a/— T • x 2 , 

b = b 1 + V^l-b 2 , m = mi + V-L ■ m 2 , y = y 1 + ■ y 2 , 

whose complex conjugates are a,n,x,b,m,y, respectively. These correspond to the vari- 
ables used in the formulation in [33]. The four quadratic polynomials of f-p are 

fi = ni - aiXi - a 2 x 2 , f 2 = n 2 + a x x 2 - a 2 xi , 
/ 3 = mi - biyi - b 2 y 2 , / 4 = m 2 + b x y 2 - b 2 y x . 

The eight quartic polynomials depend upon the displacements from Po to the other points 
Pj. For j = 1, . . . , 8, define Qj := (Qj t i, Qj, 2 ) = Pj — Po and write each displacement Qj 
using isotropic coordinates, namely (8j,8j) where 

Sj = Qj t i + v / -l" ■ Qj,2 and 5j = Q jtl - v / — T ■ Qj, 2 - 
1, . . . , 8, the quartic polynomial f i+ j of f-p is 

U+j ■= lilj + 7i7? + Ijlj 
7j := g^rj - gjrj, 7^ := r|pj - rjpj, 7° := pjgj - p^q^ 

p^ :—n — 6jX, qj := n — SjX, r X j := Sj(a — x) + Sj(a — x) — 6j6j, 

p y j:=fn-5jy, q^:=m-5jy, rjf := Sj(b -y) + 5j(b - y) - 8j8j. 

We first certified that, for nine randomly selected points in the complex plane, the 
resulting polynomial system has at least 8652 isolated solutions. Since the displacements 
Qj define the polynomial system, we choose them to be points of Q[a/— T] 2 with each 
coordinate having unit modulus of the form 

t 2 - 1 — 2t 

where t was a quotient of two ten digit random integers. We used regeneration [2T] in 
Bertini [1] to compute 8652 points that were heuristically within 1CT 100 of an isolated 
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solution for f-p. Then, alphaCertified produced a soft certificate using 256-bit precision 
that these 8652 points are approximate solutions to f-p with distinct associated solutions. 

We next certified the number of real solutions for a specific set of nine real points, 
namely Problem 3 of [33]. The nine real points are listed in Table 2 of [33], which, for 
convenience, we list the values of Sj in Tabled) Since the points are real, Sj is the conjugate 

Table 1. Values of 5, for Problem 3 of [331 



3 


Si 


1 


0.27 + 0.1^^1 


2 


0.55 + 0.7 v /r T 


3 


0.95 + ^^1 


4 


1.15 + 1.3^-L 


5 


0.85 + 1.48^^1 


6 


0.45 + 1.4^-L 


7 


-0.05 + ^^1 


8 


-0.23 + 0.4^^1 



of Sj. We used parameter continuation in Bertini to solve the resulting polynomial system 
starting from the 8652 solutions to the polynomial system solved in the first test. This 
generated a list of 8652 points which alphaCertified soft certified using 256-bit precision 
to be approximate solutions that have distinct associated solutions of which 384 are real. 
In particular, this confirms the results reported in Table 3 of [33] for Problem 3, namely, 
that 64 = 384/6 of the 1442 mechanisms are real. 

Figure [1] shows three of the 64 real mechanisms that solve this synthesis problem, to- 
gether with their workspace curves. The first has two assembly modes with the workspace 



Figure 1. Three solutions. 

curve of one mode a simple closed curve that contains the nine target points. This mech- 
anism is the only viable mechanism among the 64 real mechanisms. The second has only 
one assembly mode, but its workspace curve is convoluted and does not meet the target 
points in a useful order. The third has two assembly modes, and each only reaches a 
proper subset of the target points. 



ALPHACERTIFIED: CERTIFYING SOLUTIONS TO POLYNOMIAL SYSTEMS 15 

5.3. Lines, points, and conies. We consider geometric problems of plane conies in C 3 
that meet k points and 8 — 2k lines for k = 0, . . . , 4. When the points and lines are general, 
the numbers of plane conies are known and presented in Table [2j 

Table 2. Numbers of plane conies 



k 


4 


3 


2 


1 





Number of conies 





1 


4 


18 


92 



This problem is from a class of problems in enumerative geometry — counting rational 
curves — that has been of great interest in recent years [15] . For problems of enumerating 
rational curves of degree d in the plane that interpolate 3d— 1 real points, Welschinger [43] 
defined an invariant Wd which is a lower bound on the number of real rational curves, 
and work of Mikhalkin [31] and of Itenberg, Kharlamov, and Shustin showed that Wd is 
positive and eventually found a formula for it [23] . 

We used alphaCertified to investigate the possible numbers of real solutions to these 
problems of conies when their input data (points and lines) are real. Of particular interest 
is the minimum number of solutions that are real. Our experimental data suggests that 
when k — 1 at least two of the solutions will be real, and it shows that for k = 0,2, it is 
possible to have no real solutions. 

This experiment computed random instances of the problem. The coordinates of points 
were taken to be the quotient of two random ten digit integers, and the real lines were 
taken to be lines through two such points. The resulting polynomial system was square. 
Each real instance was solved by Bertini [4] using a straight line parameter homotopy 
starting with a fixed random complex instance (see [3H| for more details). This gave points 
that were heuristically within 10 -75 of each isolated solution. Then alphaCertified used 
256-bit precision to softly certify that the points computed by Bertini were approximate 
solutions corresponding to distinct solutions, and to count the number of real solutions. 
Since enumerative geometry provides the generic root count, this yields a post-processing 
certificate that Bertini has indeed computed an approximate solution corresponding to 
each solution to the polynomial system. 

In every instance that Bertini successfully tracked every path, the heuristic results of 
Bertini matched the certified results of alphaCertified. Out of the over 1,450,000,000 
paths tracked, 76 paths were truncated by Bertini due to a fail-safe measure. Thirty- 
two paths were truncated since they needed more than the fail-safe limit of 10,000 steps 
along the path. Each of these paths were successfully tracked when the limit was raised 
to 25,000 steps. Forty-four paths were truncated since the adaptive precision tracking 
algorithm [5], El [3] requested to use more than the fail-safe limit of 1024-bit precision. 
Each of these paths were successfully tracked when the fail-safe limit was raised to 1284- 
bit precision. 

The first interesting case is when k = 2 and there are four conies meeting two points and 
four lines. We solved 500 random real instances using the Brazos cluster. Each instance 
took an average of 0.7 seconds for Bertini to solve and 0.1 seconds for alphaCertified to 
certify the results. We found that there can be 0, 2, or 4 real solutions. Table [3] presents 
the frequency distribution of these 500 instances for this case. 
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Table 3. Frequency distribution for conies through two points and four lines 



# real 





2 


4 


total 


frequency 


12 


221 


267 


500 



When k — 1, there are 18 conies meeting a point and six lines in C 3 . We solved 
1,000,000 random real instances using the Brazos cluster. Each instance took an average 
of 1.6 seconds for Bertini to solve and an average of 0.1 seconds for alphaCertified to 
certify the results. Every real instance that we computed had at least 2 real solutions. 
Table H] presents the frequency distribution of these 1,000,000 instances for this case. 



Table 4. Frequency distribution for conies through a point and six lines 



# real 





2 


4 


6 


8 


10 


12 


14 


16 


18 


total 


frequency 





3281 


21984 


88813 


193612 


261733 


226383 


137074 


53482 


13638 


1000000 



To compare the performance of alphaCertified to symbolic methods, we computed 
40,000 instances of the conic problem with k = 1 using Singular [IT] to compute an 
eliminant that satisfies the Shape Lemma [7] and Maple to count the number of real 
roots of the eliminant, which is a standard symbolic method to determine the number of 
real solutions to a zero-dimensional system of polynomial equations. The coordinates of 
points were taken to be rational numbers p/q where p, q were integers with \p\ < 4000 and 
< q < 1000. Each computation took approximately 661 seconds on a single node of a 
server with four six-core AMD Opteron 8435 processors and 64 GB of memory. Table [5] 
presents the frequency distribution of these 40,000 instances for this case. 

Table 5. Frequency distribution for conies through a point and six lines 



# real 





2 


4 


6 


8 


10 


12 


14 


16 


18 


total 


frequency 





146 


892 


3558 


7739 


10575 


8965 


5488 


2089 


548 


40000 



Finally, when k = 0, there are 92 plane conies meeting eight general lines in C 3 . We 
solved 15,662,000 random real instances using the Brazos cluster. On average, each in- 
stance took 8.8 seconds for Bertini to solve and 0.7 seconds for alphaCertified to certify 
the results. Table |6] presents the frequency distribution of these instances. 

5.4. A Schubert problem. Our last example concerns a problem in the Schubert calcu- 
lus of enumerative geometry, which is a rich class of geometric problems involving linear 
subspaces of a vector space. Many problems in the Schubert calculus are naturally formu- 
lated as overdetermined polynomial systems. We investigate one such problem that can 
also be formulated as a square polynomial system using the approach of [2]. In particu- 
lar, we demonstrate alphaCertified's algorithms for overdetermined systems as well as 
investigate a conjecture on the reality of its solutions. 

This problem involves four-dimensional linear subspaces (four-planes) H of C 8 that 
have a non-trivial intersection with each of eight general three-planes K , . . . ,K-?. The 
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# real 





2 


4 


6 


8 


10 


12 


14 


frequency 


1 


8 


26 


65 


466 


1548 


4765 


11928 


# real 


16 


18 


20 


22 


24 


26 


28 


30 


frequency 


26439 


52875 


98129 


167932 


270267 


404918 


569891 


756527 


# real 


32 


34 


36 


38 


40 


42 


44 


46 


frequency 


942674 


1114033 


1246533 


1332289 


1355320 


1319699 


1226667 


1091019 


# real 


48 


50 


52 


54 


56 


58 


60 


62 


frequency 


932838 


762463 


596174 


449021 


323927 


223455 


149629 


95740 


# real 


64 


66 


68 


70 


72 


74 


76 


78 


frequency 


59141 


34834 


19516 


10672 


5671 


2744 


1290 


530 


# real 


80 


82 


84 


86 


88 


90 


92 


total 


frequency 


204 


90 


26 


11 


3 


2 





15662000 



Schubert calculus predicts 126 such four-planes. To formulate this Schubert problem, 
consider H to be the column space of a 8 x 4 matrix in block form 

h 
X 



H 



where I± is the 4x4 identity matrix and X is a 4 x 4 matrix of indeterminates. Represent 
a three-plane K as the column space of a 8 x 3 matrix of constants. Then the condition 
that H meets K non-trivially is equivalent to the vanishing of the determinants of the 
eight 7x7 square submatrices of the 8x7 matrix 

(9) A — [H K] . 

In this standard formulation, the Schubert problem is a system of 64 equations in 16 
indeterminates. Using a total degree homotopy to solve this would follow 4 16 paths. 
There is a second formulation which we used. Write K in block form, 



K 



1C 2 



where /Q and /C2 are 4x3 matrices. A linear dependency among the columns of A 
is given by vectors v G C 4 and w G C 3 such that Hv + Kw = 0. Applying this to the 
different blocks of H and K gives 



I4V + KL\W = and Xv + K, 2 w 







which is equivalent to Aw = 0, where A := IC2 — XK,\. Thus H meets K non-trivially 
if and only if each 3x3 minor of A vanishes. This gives a system Fq(x) of 32 cubic 
polynomials in 16 indeterminates, which is more compact than the original formulation. 

We certified solutions to this overdetermined polynomial system Fq- We randomized 
Fq to maintain the structure of the equations as follows. For each % = 0, . . . , 7 and 
j = 1,2,3,4, let fij be the determinant of the submatrix created by removing the j th 
row of the matrix A4 corresponding to the ith three-plane. Then, for each j, we take 
four random linear combinations of the polynomials /nj, f\j, . . . , fyj. This preserves the 
multilinear structure of the equations in the four variable groups corresponding to the 
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columns of X. Solving this system using regeneration [21] finds 22,254 solutions of which 
alphaCertified soft certified using 256-bit precision that 126 of these are approximate 
solutions to two different random square subsystems of Fo with associated solutions within 
a distance of 5 = 10~ 10 of each other. The same result was also obtained using 6 = 10~ 5 . 
Thus, alphaCertified provided a soft certificate based on the heuristic algorithm for 
overdetermined systems that we found all 126 solutions to the Schubert problem. 

This Schubert problem has an equivalent formulation as a square system. The columns 
of A are linearly dependent if and only if there exists ^ v E C 3 such that Av = 0. For 
generic ati, a 2 E C, this occurs if and only if there exists yi, y 2 E C such that 



A ■ 



Hi 

a 1 y 1 + a 2 y 2 + 1 



This yields a system of 32 polynomials in 32 indeterminates, say F$(x, y^ . . . , y^ >). This 
polynomial system consists of four bilinear polynomials in x and for each i = 0, . . . , 7. 
Since consists of two indeterminates, namely y^ and y 2 \ a 9- homogeneous homotopy 
used to solve Fs would follow ( 4 ) = 6 8 paths. As described in [2], we are interested in the 
components of V(Fs) having fibers with generic dimension zero. For generic Kq, . . . , Kj, 
since V(Fs) is zero-dimensional, V(i*b) and V(Fs) both consist of 126 isolated points and 
V(F S ) naturally projects onto V(F G ). 

We also investigated the number of real solutions when the three planes Ki are as 
follows. For t E M, let ^(t) = (1, t, t 2 , . . . , t 7 ) G K 8 be a point on the moment curve. Select 
24 rational numbers t 1; . . . , t 2 4 and for % — 0, . . . , 7, let Ki be the span of the three linearly 
independent vectors 7(t 3 j + i), 7(^+2), and 7(^+3). When t x < t 2 < • • ■ < t 2 4, the Secant 
Conjecture [TH] posits that all 126 solutions will be real, but if the points are not in this 
or some equivalent order, then other numbers of real solutions are possible. 

Since Kq, . . . , K-j are real, if the constants on are real, then the real points of V(Fo) 
correspond to the real points of V(Fs). We solved 25000 real instances using random 
numbers in the interval [—1,1] and softly certified that each had 126 real solutions, using 
256-bit precision. 

Lastly, we investigated the number of real solutions when the three planes Ki are as 
follows. For z = 0, . . . , 7, let tj E C be generic under the condition that 2k are complex 
conjugate pairs and 8 — 2k are real, where < k < 4. Define Ki = T(ti) where 



T(t) 



- 1 





" 


t 


1 





t 2 


2t 


1 


t 3 


3t 2 


3t 


t 4 


4t 3 


6t 2 


t 5 


5t 4 


10t 3 


t 6 


6t 5 


15t 4 



Then Ki is the three-plane osculating the moment curve at the point 7(tj). When k = 0, 
that is, when each ti is real, this is the Shapiro Conjecture (MTV Theorem) [3U El] 
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and all 126 solutions are real. We tested 1000 such instances and for each, alphaCer- 
tified correctly identified all 126 solutions to be real. Our primary interest was when 
k > 0, for we wanted to test the hypothesis that there would be a lower bound to 
the number of real solutions if the set of osculating three-planes were real (that is, if 
{Kq, . . . , K7} = {Kq, . . . , Kj}). This is what we found, as can be seen in the partial fre- 
quency table we give below. (To better show the lower bounds, we omit writing in the 
cells with no observed instances.) This enumeration of real solutions was softly certified 
using 256-bit precision. 

Table 7. Frequency distribution for the Schubert problem 





# real 






k 





2 


4 


6 


8 


10 


12 




18 


20 


22 




124 


126 


total 































1000 


1000 


1 








6 


6 


10 


88 




554 


1888 


1832 




69 


2021 


42000 


2 












2614 


3771 




3285 


1579 


1378 




1 


38 


24000 


3 












8896 


4479 




1079 


721 


2586 








23500 


4 






















19134 






1 


22500 



This computation was part of a larger test of hypothesized lower bounds [20] . 

6. Conclusion 

Smale's a-theory provides a way to certify solutions to polynomial systems, determine 
if two points correspond to distinct solutions, and determine if the corresponding solu- 
tion is real. Using either exact rational or arbitrary precision floating point arithmetic, 
alphaCertified is a program which implements these a-theoretical methods. 

We have also produced a Maple interface to alphaCertified to facilitate the construc- 
tion of the input files needed. 
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